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ABSTRACT 


Two-dimensional  numerical  solutions  of  three  problems  of  stress 
wave  Interactions  with  oblique  cracks  in  granite  have  been  carried  out  using 
the  SHBP  code,  a  finite-difference  Lagrangian  program  that  incorporates  a 
comprehensive  hyclrodynainic-elastic-plastic  behavioral  model.  The  problems 
involve  both  Infinite  and  finite-length  crack  surfaces.  To  enable  the  study 
of  propagation  of  cracks  und.r  stress  wave  loading,  a  dynamic  Griffith  criterion 
has  been  formulated  for  incorporation  into  the  code.  Comparisons  of  the 
numerical  results  with  analytic  solutions  for  model  problems  are  being  obtained 
to  verify  the  code.  Included  in  these  correlations  is  the  problem  of  an  accel¬ 
erating  crack  for  tiie  case  of  anti-plane  shear. 
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1. 


INTRODUCTION 


The  objective  of  this  program  is  to  apply  two-dimensional  numerical 
techniques  to  the  solution  of  problems  of  stress  wave  interactions  with  indi¬ 
vidual,  oblique  cracks  in  rock  media.  Computer  codes  suitable  for  solving 
a  wide  range  of  fluid  and  solid  mechanics  problems  have  been  c  vailable  or 
under  development  for  several  years  and  efforts  to  extend  these  techniques  to 
the  quantitative  analysis  of  wave  interactions  with  cracks  and  of  crack  motion 
are  now  feasible.  Shock  Hydrodynamics  two-aimensional  SHEP  (Shock 
Hydrodynamic  _Elastic_Plastic)  code  is  being  utilized  for  this  purpose  in  this 
program. 

This  report  describes  the  work  conducted  during  the  first  six  months 
of  the  program. 

2.  SUMMARY 

2.1  PROBLEM  AREA 

Excavation  processes  in  rock  media  typically  involve  the  action  of 
strong  dynamic  stresses  Introduced  either  from  explosive,  mechanical,  or 
other  impulsive  loading  sources.  The  propagation  of  stress  waves  in  homo¬ 
geneous,  isotropic  media  is  reasonably  well  understood.  In-situ  rock  media, 
however,  typically  contain  large  scale  discontinuities  in  the  form  of  cracks, 
joints,  and  faults.  Interactions  of  stress  waves  with  these  discontinuities 
can  cause  slippage  along  '.he  cracks,  or  extension  (prop?  go  .ion)  of  the  cracks, 
or  separation.  The  interactions  can  also  alter  the  ch  iractnristics  of  the  stress 
wave  transmitted  across  the  discontinuity. 

This  study  is  concerned  with  an  analysis  of  the  detailed  mechanisms 
involved  when  strong  stress  waves  interact  with  the  crack  surfaces  in  jointed 
block  media.  Of  particular  interest  are  cracks  which  are  obliquely  oriented 
relative  to  the  wave  front.  The  understanding  thus  obtained  can  contribute 
to  the  advancement  of  knowledge  of  excavation  processes  in  various  media, 
and  how  to  control  and/or  improve  such  processes. 
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The  technical  approach  that  is  being  used  in  studying  the  details  of 
stress  wave-crack  interactions  is  based  on  two-dimensional  numerical  analyses 
of  the  dynamic  phenomena  occurring  under  various  conditions  of  stress  wave 
profile,  crack  orientation  relative  to  the  wave  front,  and  degree  of  locking 
(or  lubrication)  initially  fc  und  across  the  crack.  The  computer  program  being 
used  to  obtain  the  numerical  solutions  is  the  SHEP  code.  SHEP  is  a  finite- 
difference  Lagrangian  program  employing  a  comprehensive  hydrodynamic- 
elastic-plastic  behavioral  model.  SHEP  has  been  under  intensive  use  and 
development  for  the  past  six  years  and  has  been  applied  to  a  broad  spectrum 
of  wave  propagation  problems. 

A  major  difficulty  in  examining  wave  interactions  with  discontinuities 
such  as  cracks  or  fracture  surfaces  arises  due  to  the  constraint  of  the  contin¬ 
uum  model  which  is  normally  assumed  in  numerical  analyses  of  wave  propaga¬ 
tion.  Special  routines  in  the  SHEP  code  alleviate  this  difficulty  by  permitting 
crack  surfaces  to  be  explicitly  defined  in  the  computational  grid.  Thus  the 
grid  is  not  coupled  across  the  crack,  and  slippage  and/or  separation  can 
occur. 


2.2  PLAN  OF  RESEARCH 

The  work  under  the  current  contract  is  divided  into  the  following  two 
tasks,  corresponding  to  analyses  of  stress  wave  interactions  with 

a)  single,  infinite  crunks  and 

b)  single,  finite-length  cracks. 

2,2.1  Task  1  -  Interaction  of  Stress  Waves  with  Single, 

Infinite  Cracks 


Tils  task  is  concerned  with  the  analysis  pertaining  to  cracks  which 
are  infinit  in  extent,  or  which  intersect  with  tie  ground  surface.  The  work 
initially  consists  of  the  selection  and  determination  of  the  problem  specifi¬ 
cations  (such  as  media  properties,  loading  wave  characteristics,  and  crack 
orientation  and  condition).  SHEP  code  solutions  of  the  problems  defined  are 
then  sot  up  and  run.  These  solutions  provide  complete  quantitative  data  for 
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all  the  state  and  motion  variables  of  interest  throughout  the  computing  field 
at  regular  intervals  of  time.  In  addition,  spatial  plots  of  the  principal  stress, 
particle  velocity,  and  displacement  fields  are  obtained  at  selected  times 
during  the  event.  Following  completion  of  the  solutions,  analysis  and  inter¬ 
pretation  of  the  results  are  performed  *o  assess  the  important  mechanisms  and 
factors  influencing  the  behavior  of  the  crack  and  the  characteristics  of  the 
transmitted  stress  wave. 

2.2.2  Task  2  -  Interaction  of  Stress  Waves  with  Single , 

Finite  Cracks 


This  task  is  similar  to  Task  1,  except  that  the  class  of  cracks 
involved  are  those  of  finite  length,  so  that  the  interaction  of  a  stress  wave 
with  a  crack  tip  can  be  examined.  Particular  attention  is  being  given  to  the 
stress  magnitudes  developed  near  the  tip,  and  the  effect  of  the  tip  on  the 
transmitted  wave  system.  The  problems  are  selected  so  that  comparisons 
between  the  cases  of  infinite  (Task  1)  and  finite  (Task  2)  cracks  can  be  made. 

Also  included  in  this  task  is  the  development  of  a  dynamic  Griffith 
criterion  to  be  incorporated  into  the  code  for  the  study  of  crack  propagation 
under  stress  wave  loading. 

To  verify  the  suitability  and  accuracy  of  the  code,  and  modifica¬ 
tions  made  thereto  during  the  course  of  the  program,  for  analyses  of  wave/ 
crack  interactions,  analytic  solutions  of  selected  test  problems  are  being 
obtained  and  used  to  check-out  corresponding  numerical  solutions  obtained 
with  the  code. 

2.3  MAJOR  ACCOMPLISHMENTS 

Primary  program  accomplishments  during  this  period  included 

a)  The  completion  of  three  SHEP  code  solutions  of  sL.ess 
wave/crack  interaction  problemj, 

b)  the  initial  formulation  of  a  dynamic  Griffith  criterion 
for  crack  propagation,  suitable  for  incorporation  Into 
the  code,  and 
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c)  the  initiation  of  a  series  of  test  cases  to  verify  the 
numerical  results  through  comparison  with  analytic- 
solutions  that  are  being  obtained. 

2.3.1  Numerical  Solutions 

Three  problems  of  wave  interactions  with  cracks  were  selected  for 
analysis  by  means  of  code  solution,  as  depicted  in  Figure  1.  These  problems 
were  chosen  to  demonstra:e  the  utility  of  numerical  techniques  for  obtaining 
detailed  information  on  the  response  of  cracked  media  subjected  to  impulsive 
loads,  in  general,  and,  in  particular,  for  assessments  of  the  wave  interactions 
in  the  vicinity  of  a  crack.  It  is  noted  that  the  explicit  definition  of  a  crack 
surface  such  as  specified  in  these  problems  is  not  generally  amenable  to 
treatment  through  conventional  code  techniques,  which  normally  assume  a 
continuous  material  model. 

Tne  rock  medium  selected  for  these  problems  was  granite.  The 
material  properties  assumed  for  the  granite  are  described  in  Section  3.2  of 
this  report.  The  problems  were  run  in  plane  geometry,  assuming  plane  strain; 
the  variables  are  thus  independent  of  the  z-coo; dinate  (perpendicular  to  the 
cross-section  shown). 

The  first  problem  selected  for  analysis  (Case  1)  consisted  of  inter¬ 
actions  of  a  stress  wave  with  a  crack  oriented  at  a  30°  angle  with  the  wave 
front.  The  stress  wave  was  generated  by  uniformly  loading  the  left  face  of 
the  granite  block  with  a  pressure  pulse.  A  triangular  pressure  pulse  of  5 
kllobars  peak  magnitude  and  0.2  millisecond  duration  was  used  for  this 
problem,  as  sketched  below: 
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CASE  2.  INTERACTION  OF  STRESS  WAVE  WITH  SINGLE,  FINITE- LENGTH  CRACK 


CASE  3.  INTERACTION  OF  STRESS  WAVE  WITH  SINGLE,  FINITE -LENGTH 
CRACK,  WITH  CRACK  GROWTH 


Crack  Tip 


Figure  1.  Specifications  of  Problems  for  SHEPCode  Solutions. 
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The  crack  in  this  problem  was  characterized  by  a  free  slip  condition 
and  zero  width.  Opening  of  the  crack  was  allowed  to  occur  if  stress  compon¬ 
ents  normal  to  the  crack  went  into  tension. 

Tne  second  problem  considered  (Case  2)  involved  the  interaction  of 
a  stress  wave  with  a  finite-length  crack.  The  angle  of  orientation  of  the 
crack  and  the  applied  pressure  loading  were  the  same  as  in  the  first  problem. 
The  crack  extended  from  the  loading  surface  (lower  left)  to  the  crack  tip, 
situateci  on  the  horizontal  mid-plane  of  the  block. 

Case  3  was  the  same  as  the  second  problem,  except  that  crack 
growth  was  permitted.  This  case  demonstrates  the  provisions  in  the  code 
which  can  be  used  to  model  crack  propagation.  In  this  case,  dynamic  decoup¬ 
ling  of  lattice  points  in  the  computing  mesh  occurs  when  a  specified  criterion 
is  satisfied. 

SHEP  code  solutions  of  these  three  problems  were  successfully 
completed.  Plots  depicting  the  particle  velocity  field  and  the  principal  stress 
field  occurring  in  the  test  block  were  obtained  for  several  times  during  the 
interactions.  In  addition,  time  histories  of  pertinent  parameters  at  several 
stations  in  the  field  were  recorded.  These  results  are  discussed  in  detail 
in  Section  3  of  the  report.  Some  representative  results  of  these  code  solu¬ 
tions  are  shown  here,  in  Figures  2  to  7. 

For  Case  1,  the  principal  stress  field,  for  a  time  of  .3  msec,  and 
the  particle  velocity  field,  for  a  time  of  .5  msec,  are  shown  in  Figures  2  and 
3.  As  the  wave  encounters  the  crack,  the  principal  stress  vectors  may  be 
seen  (Figure  2)  to  rotate  into  a  direction  transverse  to  the  crack  surface, 
reflecting  the  fact  that  the  crack  surface  can  not  bear  shear  stress.  This 
interaction  produces  a  dilatational  wave  and  trailing  shear  wave  which  propa¬ 
gate  act oss  the  block,  as  indicated  in  the  velocity  field  plot  (Figure  3).  The 
peak  stress  in  the  transmitted  wave  was  reduced  by  about  25%  from  that  in 
the  incident  wave. 

An  example  of  the  results  of  the  Case  2  solution  is  shown  in 
Figure  4,  which  depicts  the  particle  velocity  field  occurring  at  a  time  of 
.6  msec.  In  this  interaction,  the  wave  system  is  divided  approximately  in 
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Figure  3.  Particle  Velocity  Field,  Case  1 
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Figure  5.  Particle  Velocity  Field,  Case 
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Figure  6.  Particle  Velocity  Field,  Case 
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half,  the  part  above  the  crack  tip  appearing  as  a  simple  plane  wave,  and 
that  below  as  a  dilatational  wave  and  trailing  sheoi  wave,  as  in  the  previous 
solution.  Starting  from  the  crack  tip,  a  disturbance  propagates  into  the 
plane  wave  region  above  and  the  "cracked"  region  below,  altering  both  flow 
fields  and  creating  an  expanding  region  of  transition  between  them. 

In  Case  3,  where  crack  growth  was  allowed,  the  criterion  used  for 
decoupling  of  points  beyond  the  crack  tip  was  that  each  of  the  cells  surround¬ 
ing  a  lattice  point  must  have  failed,  i.e. ,  at  some  time  reached  a  state  on 
the  granite  failure  surface.  This  was  a  conservative  criterion,  since  the 
failure  surface  generally  represents  states  where  virtually  complete  fracture 
occurs.  Use  of  more  sensitive  criteria  can  be  employed,  and  development 
of  a  dynamic  Griffith  criterion  is  currently  under  development,  as  described 
below.  The  crack  growth  occurring  in  this  solution  is  indicated  in  Figures  5, 
6,  and  7,  which  are  plots  of  the  particle  velocity  field  for  times  of  .5,  .6, 
and  .7  msec.  The  extent  of  the  crack  in  each  plot  is  indicated  by  the  circled 
lattice  points. 

2.3.2  Dynamic  Griffith  Criterion 

In  connection  with  efforts  being  made  under  the  program  to  enable  the 
study  of  propagation  of  cracks  under  stress  wave  loading,  the  incorporation  of 
equations  into  the  SHEP  code  which  govern  the  rate  of  propagation  of  a  brittle 
crack  surface  in  an  elastic  material  is  currently  under  development.  These 
equations  are  known  as  the  Griffith  criterion  and  they  provide  a  relation 
between  power  input  to  the  body  and  the  rate  of  uptake  of  this  power  by  strain 
energy,  kinetic  energy,  and  new  surface  energy.  The  concept  of  surface 
energy  is  the  feature  that  was  introduced  by  Griffith  in  the  early  1900's,  and 
it  requires  the  determination  of  an  additional  material  parameter,  namely  the 
surface  energy  per  unit  area.  It  is  planned  to  extract  the  value  of  this  para¬ 
meter  for  granite  from  available  fracture  data. 
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2.3.3  Analytical  Comparison  Problems 

To  verify  the  SHEP  code  solutions  and  the  formulation  changes  being 
made,  comparisons  of  numerical  results  with  analytical  solutions  of  model 
problems  are  being  made. 

As  part  of  this  effort,  the  capability  has  been  added  to  the  code  for 
the  treatment  of  anti-plane  shear,  or  out-of-plane  displacements,  with  the 
restriction  that  the  motions  are  independent  of  he  z-coordinate ,  so  as  to 
retain  the  two-dimensional  character  of  the  code.  This  was  done  primarily 
since  the  only  elasto-dynamic  solutions  currently  available  for  an  accelera¬ 
ting  crack  are  those  for  the  case  of  anti-plane  shear,  although  it  also  repre¬ 
sents  a  useful  tool  in  numerical  analysis  which  has  heretofore  been  unavailable. 

A  model  problem  of  simple,  shear  motion  of  a  slab  has  been  solved 
with  the  modified  code.  The  results  of  the  code  solution  showed  excellent 
agreement  with  the  analytical  solution  for  this  case.  The  next  case  in  this 
series,  currently  in  work,  is  a  problem  involving  the  interaction  of  an  anti¬ 
plane  shear  wave  with  a  stationary  crack. 

3.  NUMERICAL  SOLUTIONS 


As  noted  above,  numerical  solutions  of  three  problems  involving  the 
interaction  of  stress  waves  with  cracks  were  performed.  The  specifications 
of  these  problems  were  given  in  Section  2.3.1. 

3. 1  COMPUTATIONAL  METHOD 

3.1.1  Physical  Model 

The  computer  program  being  used  in  this  study  is  the  two-dimensional 

SHEP  code,  which  solves  the  equations  of  motions  for  elastic-plastic  bodies 

by  means  of  a  finite-difference  Lagrangian-cell  technique.  SHEP  has  been 

under  intensive  use  and  development  for  the  past  six  years  and  has  been 

previously  documented*  and  distributed  to  interested  parties.  The  mathema- 

2 

tical  formulation  is  basically  the  same  as  that  described  by  Wilkins  .  To 
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delineate  the  boundary  between  elastic  and  plastic  deformations,  various  yield 
criteria  may  be  used,  such  as  von  Mises,  Mohr -Coulomb,  or  arbitrary  functions. 
Within  the  chosen  yield  surface,  the  deformations  are  considered  to  bo  elastic, 
i.e. ,  when 

/ 3  V2  <  Y  (1) 

where  JJ,  is  the  second  invariant  of  the  deviatoric  stress  tensor  and  Y  is  the 
yield  strength.  Excursions  on  the  yield  surface  arc  made  in  accordance  with 
the  Prandtl-Reuss  flow  rule. 

To  model  a  crack  surface,  SHEP  contains  provisions  for  inserting 
surfaces  of  discontinuity,  which  consist  of  grid  lines  having  a  dual  set  of 
lattice  points.  These  surfaces  are  discussed  in  the  following  section. 

3.1.2  Surfaces  of  Discontinuity 

In  a  normal  Lagrangian  computational  grid,  material  elements  on 
either  side  of  an  interface  at  any  point  are.  coupled  to  each  other  for  the  entire 
problem;  they  are,  i.e. ,  locked  or  welded  together  along  the  line  segment 
connecting  any  two  lattice  points  along  the  interface.  At  any  interface,  which 
may  represent  a  crack  within  a  material  or  the  boundary  between  two  different 
materials,  there  are,  however,  in  general,  special  boundary  conditions  which 
apply,  and,  in  addition,  there  is  the  possibility  of  forces  which  may  be  set 
up  that  tend  to  cause  the  materials  to  slip  past  each  other  or  to  separate.  A 
gas  flowing  past  a  metal  surface  is  an  example  of  such  a  case.  The  onset  of 
material  fracture  during  a  problem  also  gives  rise  to  the  requirement  for  treat¬ 
ing  the  decoupling  or  uncoupling  of  elements  which  are,  in  this  case,  within 
an  originally  competent  ma.erial.  For  application  to  problems  in  fracture 
mechanics,  such  as  in  this  program,  the  latter  requirement  is  particularly 
important. 

A  formulation  of  sliding  interfaces  for  Lagrangian  codes,  as  reported 
2 

by  Wilkins  ,  provided  a  capability  for  the  numerical  treatment  of  problems 
involving  sliding  of  tw'o  materials  along  an  interface.  This  formulation  served 
as  the  basis  for  development  of  the  surface  of  discontinuity  capability  currently 
available  in  the  SHEP  code. 


15 


The  basic  features  of  the  surface  of  discontinuity  formulation  are 
illustrated  in  Figure  8.  The  grid  line  corresponding  to  the  surface  of  discon¬ 
tinuity  is  known  in  common  parlance  as  a  slide  line.  At  the  start  of  a  problem, 
the  lattice  points  along  the  slide  line  may  be  individually  designated  as 
decoupled  points,  corresponding  to  their  lying  on  an  interface,  or  as  eoupled 
points,  in  which  case  their  behavior  is  the  same  as  in  an  ordinary  mesh.  For 
deco.’p’ed  points,  special  sets  of  governing  equations  are  used  to  individually 
determine  the  motion  of  the  point  pairs,  to  reflect  the  fact  that  there  is  an 
interface,  such  that,  e.g.,  shear  stress  cannot  be  supported.  If  forces  are 
present  which  tend  to  cause  slippage,  the  decoupled  points  will  thus  disen¬ 
gage  end  move  separately  along  the  slide  line. 

Additionally,  the  development  of  tensile  stresses  normal  to  an  inter¬ 
face  will  tend  to  cause  material  separation  and  formation  of  voids.  Provisions 
have  been  made  in  the  code  to  treat  this  phenomenon,  also. 

The  void  opening  test  is  made  by  computing  the  stress  normal  to  the 
interface  at  a  decoupled  slide  p-  .t  and  comparing  this  value  with  a  selected 
critical  value  of  stress  required  for  uncoupling.  If  the  computed  stress  is 
greater  than  the  critical  value  (in  tension),  then  that  point  is  designated  as 
a  free  surface  point.  The  newly  formed  free  point  is  then  moved  in  accordance 
with  the  regular  equations  of  motion  for  a  point  on  a  free  surface. 

For  the  problems  performed  in  this  study,  the  critical  value  for 
uncoupling  was  set  to  zero,  such  that  any  tensile  stress  would  tend  to  cause 
separation.  In  other  applications ,  e.g.,  for  the  interfaces  in  laminated 
materials,  this  value  could  be  set  equal  to  the  bond  strength. 

Void  closure  may  also  occur  and  is  treated  in  the  code  by  appro¬ 
priate  tests  to  determine  if  the  materials  have  come  in  contact.  If  so,  the 
equations  of  the  surface  of  discontinuity  are  restored,  and  the  materials  may 
subsequently  slip  or  re-open,  as  before. 

Lattice  points  along  the  slide  line  which  are  initially  designated  as 
coupled  points  may  dynamically  decouple,  individually,  during  the  course 
.f  a  problem,  if  a  selected  criterion  is  met.  Various  decoupling, or  fracture, 
criteria  may  be  used.  Once  the  lattice  points  are  decoupled,  the  equations 
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a.  Normal  Coupled 
Lagrangian  Grid 


b.  Decoupled  Grid  With  Relative 
Motion  (Slippage)  Along  Material 
Interface,  Fracture  Surface,  or 
Plane  of  Weakness 


"Crack  Propagation"  With 
Subsequent  Slippage  Along 
Fracture  Surface 


d.  "Crack  Propagation"  With 
Subsequent  Separation  and 
Rejoining  Along  Fracture  Surface 


Figure  8.  Schematic  Illustrations  of  the  Treatment  of  Slippage, 
Fracture,  Void  Opening,  and  Void  Closing  Along 
Surfaces  of  Discontinuity  in  the  SHEPCode. 
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of  the  surface  of  iiscontinu'ty  are  Invoked.  A  mechanism  is  thus  provided 
which  can  be  used  to  mode,  crack  propagation  within  a  material. 

Additional  information  on  the  code  operations  and  the  mathematical 

2-4 

brills  of  sliding  interfaces  has  been  previously  reported. 

3 . 2  MATERIAL  PROPERTIES 

The  rock  medium  in  these  problems  was  granite.  The  properties 
selected  for  the  granite  were: 

3 

Density:  Po  =2.69  gm/cm 

Dilatational  Velocity:  v^  =  .  579  cm/psec 

o 

Shear  Velocity:  v  =  .  330  cm/psec 

so  ' 

These  values  imply  the  following  other  properties: 

Bulk  Modulus:  KQ  =  .512  Mb 

Shear  Modulus:  GQ  =  .293  Mb 

Poisson's  Ratio:  v  =  .26 

o 

The  subscript  o  in  the  above  indicates  that  these  are  normal,  pre-shocked 
values.  These  values  were  selected  from  previous  studies  involving  granite 
media  (References  5  and  6). 

The  equation  of  state  of  granite,  suitable  for  the  low-pressure 
regime  applicable  in  these  problems,  was  formulated  as  follows. 

P  =  Ap  +  Bp2  6  +  Gpe  P  <  .04  Mb  (2) 

where 

6  =  1  for  p  >  0 
6  =  0  for  p  s  0 

The  symbols  are  defined  as 

e  =  specific  internal  energy 
P  =  pre  s  s  ure 
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t)  =  p/pc  =  relative  density 

u  —  T)  -  1  -  compression 
p  «s  density 

The  values  of  the  coefficients  are: 

A  =  .512  Mb 

B  =  1.49  Mb 
G  =  2.1 

No  hydrostatic  tension  was  permitted,  i.e. ,  Pmln  =  0. 

A  Mohr-Coulomb  type  yield  model  war-  used,  i.e. , 

Y  =  .0003  +  (l-e'P//*0003)  (.00094  +  1.33  P)  (3) 

where  Y  is  the  yield  strength,  in  megabars.  The  maximum  value  permitted 
for  Y  was  10  kilobars. 

3.3  CASE  1  -  INTERACTION  OF  STRESS  WAVE  WITH  SINGLE, 

INFINITE  CRACK 

The  computational  grid  set  up  for  the  code  solution  of  Case  1  is 
shown  in  Figure  9.  This  grid  contains  2690  cells,  with  the  basic  cell  size 
set  at  10  cm  x  10  cm.  Beyond  a  central  region  of  interest  the  cell  dimensions 
geometrically  increase  in  order  to  conserve  the  total  cell  count  and  computa¬ 
tional  time.  Representative  results  of  the  code  solution,  as  depicted  by 
particle  velocity  fields  and/or  principal  stress  fields  for  times  of  .3,  .5, 
and  .92  msec,  are  shown  in  Figures  10  to  12  and  in  Figures  2  and  3  in  the 
Summary,  Section  2.3.1.  For  clarity  in  reading  these  plots,  the  field  of 
view  was  limited  to  the  central  region  of  Interest. 

For  the  stress  field  plots,  the  principal  components  of  the  stress 
tensor  for  each  cell  are  shown,  as  follows:  The  magnitude  of  the  two  princi¬ 
pal  stresses  in  the  x-y  plane  are  plotted  in  their  corresponding  principal 
directions.  The  third  principal  stress  (in  the  z  direction)  is  plotted  along 
the  line  bisecting  the  other  two  principal  directions.  Vectors  pointing  to 
the  right  are  compressive,  to  the  left,  tensile.  An  example  of  how  a  stress 
tensor  is  plotted  is  sketched  below: 
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figure  9.  Initial  Configuration  of  the  Lagrangian  Computational  Grid, 
Case  1 
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Figure  11.  particle  Velocity  Field,  Case  1 
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Figure  12.  Principal  Stress  Field,  Case  1 
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The  edits  of  the  velocity  vector  field  plot  the  direction  and  magni¬ 
tude  of  the  velocity  of  each  lattice  point  in  the  computing  grid.  The  vector 
lengths  for  both  the  principal  stress  and  velocity  fields  are  scaled  to  the 
unit  length  indicated  above  each  plot;  the  units  are  Mb/cm  and  (cm//isec)/cm, 
respectively.  Op  the  more  recently  produced  plots,  the  scale  is  also  graphi¬ 
cally  depicted  in  the  upper  right  hand  corner. 

As  the  stress  wave  interacts  with  the  crack,  the  material  on  the 
right  side  of  the  crack  is  driven  by  the  stress  component  normal  to  the  crack, 
since  shear  stresses  cannot  be  supported.  As  shown  in  Figure  10,  the  velo¬ 
city  vectors  along  the  crack  at  the  shock  front  are  thus  directed  normal  to  the 
crack,  turned  downward  30°.  Also,  as  shown  in  Figure  2,  note  that  the  prin¬ 
cipal  stress  tensors  in  the  material  along  the  right  side  of  the  crack  are  rotated 
into  a  direction  transverse  to  the  crack  surface,  again  reflecting  the  fact  that 
the  crack  surface  can  not  bear  shear  stress.  As  the  incident  wave  runs  along 
the  crack,  a  dilatational  wave  and  trailing  shear  wave  are  formed  which 
propagate  across  the  block.  The  transmitted  shock  front  remains  approxi¬ 
mately  planar  and  oriented  at  90°  to  the  x  axis.  In  the  shear  region  behind 
the  shock,  a  distinctly  downward  velocity  flow  is  evident.  This  action 
Induces  material  slippage  along  the  crack,  the  material  on  the  right  side  of 
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the  crack  moving  downward  and  to  the  left,  along  the  crack,  relative  to  the 
material  on  the  left  side.  In  addition,  there  was  a  slight  separation,  or 
opening-up,  of  the  materials  on  either  side  of  the  crack,  which  were  initially 
in  contact. 

Time  histories  of  the  material  displacement  at  the  points  indicated 
in  Figure  13  were  recorded  during  the  code  solution.  Tne  slippage  of  material 
initially  at  the  point  x  =  109  cm,  y  =  0  cm,  as  given  by  the  distance  between 
points  on  opposite  sides  of  the  crack  (points  Band  C  in  Figure  13),  is  shown 
in  Figure  14.  The  extents  of  the  downward  (y-direction)  displacements  of 
these  points  are  shown  in  the  time  histories  given  in  Figure  15.  The  down¬ 
ward  displacements  of  other  points  in  the  field,  at  x  =  55,  158,  203,  and 
253  cm  (points  A,  D,  E,  and  F)  along  the  central  horizontal  plane  (y  «*  0),  are 
shown  in  Figure  16.  The  forward  (x-direction)  displacements  of  these  points 
during  this  time  were  all~  0.25  cm.  The  perturbation  of  flew,  as  measured 
by  the  downward  thrust  of  material,  is  seen  to  diminish  as  the  distance  from 
the  crack  increases. 

Stress  (o^)  -  time  profiles  at  points  A,  D,  E,  and  F  are  shown  in 
Figure  17.  The  peak  stress  in  the  transmitted  wave  (points  D,  E,  and  F)  is 
seen  to  be  reduced  by  about  25%  from  that  in  the  incident  wave.  The  afore¬ 
mentioned  two-wave  structure  in  the  transmitted  wave  and  the  reflected  wave 
at  point  A  are  also  displayed  in  these  plots. 

3.4  CASE  2  -  INTERACTION  OF  STRESS  WAVE  WITH  SINGLE, 

FINITE-LENGTH  CRACK 

The  initial  configuration  of  the  Lagrangian  grid  set  up  for  this 
problem  is  shown  in  Figure  18.  The  crack  extends  from  the  loading  surface 
(lower  left)  to  the  crack  tip  at  x  =  200  cm,  /  =  0. 

Representative  results  of  the  SHEP  code  solution  of  this  problem, 
as  depicted  by  the  particle  velocity  fields  for  times  of  .3,  .4,  .5,  .6,  and 
1  msec,  are  shown  in  Figures  19  to  22  and  in  Figure  4  in  the  Summary, 

Section  2.3.1.  Associated  principal  stiess  fields  for  times  of  .4,  .5,  and 
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Figure  13.  Location  of  Time  History  Data  Stations 
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Figure  16.  Vertical  Displacement  at  Four  Stations  Along  the 
y  =  0  Plane,  Case  1. 
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Plane,  Case  1. 
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Initial  Configuration  of  the  Lagrangien  Computational  Grid 
Case  2 
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Figure  20.  Particle  Velocity  Field,  Case  2. 


T  ,r  K>  °  °  n 

O  1  ’  _  O  •_  °Q 

O  CJ 

°  e>  °»>  n 

s{  * 

V 

c 

51 

o  a  b0  a  b  ■  „  u„  ■„  •*»  b 

°"«.  ■  o  CO 

CD 

a  ■< 

CO 
O  * 

•  o 

1 

t 

•  • 

°  •  • .  •  * 
1  S  ,1. tVl*.! 

o  ° 

^  uo 

II 

1— 

•  •  *  t 

•  •  • 

•  •  •  • 

•  •  •  • 

1  ■  *  .  *  . '  »  .  o 

I  iiif 

»  ►  *  |i  U I  * 
ft 5  Vf1  ^  1 

o  f  ^  ™  f!u M mwfjl  $ 

BiPlJB! 

i  v,1!  'i  iv;^, 1  v  $ 

n 

*  •*  *  1 ,  »  1 
c  .  •  *  • 

o  J  •  •  . 

s  5 

**  f  #  • 

*  %  \  ’  f  ' 

£  .  •  *.  » 

e  • 

fc 

'  <V'V  -'v  ''l'"- >-L 
o /’-/V -0-V--  N 

•  •  i' •  fH*  ** 

•V  «’,> > ; 

*  • .  •  • 

•-S;  •  V  •*  o 

-  ~v.  a  *  .  - 

>•  .  -  •».  *  n 

5  .  ■  /  . 

•  .  ,  • 

r»  .  .  *  • 

•  ,  <  .  * . 

— r— .  v  s,  > 

•  .*  .  *.•.*.  • .  •  ,  •  .*  «  ' ,  '  /  ' 

•  •  •  .  *  •  #  b  #  '  .  % 

•  *  i  » .  •  •  *  .  •  .  •  ' .  ' 

.  * .  •  * .  * .  • ,  •  * .  *.-.*%  • 

•  *.  •/•/-  *.  ••  •/•  *.  *.*  •.* 

/  »  s  •_  O 

•  >  /  . 

V  *V  /  • 

/ 

»  *  ,  •  ' .  * 
i  ,  •  •  *• 

SilVVV 

?2t 

<*  *  •  * m  •  a  <*  « 

* y  yT  ° 
j  x 

0*002  0* 00T  0  0  001-  0'002 

WO  A 


34 


l 


Figure  21.  Particle  Velocity  Field,  Case  2 
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Figure  22.  Particle  Velocity  Field,  Ca 


.7  msec  are  shown  in  Figures  23  to  25.  For  clarity  in  reading  these  plots, 
the  field  of  view  is  limited  to  the  central  region  of  interest. 

The  response  of  the  granite  in  this  problem  is,  of  necessity,  similar 
to  that  in  the  infinite-crack  problem,  until  the  wave  front  reaches  the  crack 
tip.  Subsequently,  the  wave  system  is  divided  approximately  in  half,  the 
part  above  the  crack  tip  appearing  as  a  simple  plane  wave,  and  that  below 
as  a  dilational  wave  and  trailing  shear  wave,  as  in  the  previous  solution. 
Starting  from  the  crack  tip,  a  disturbance,  or  bow  wave,  propagates  into  the 
plane  wave  region  above  and  the  "cracked"  region  below,  altering  both  flow 
fields  and  creating  an  expanding  region  of  transition  between  them.  The 
diversion  of  flow  around  the  crack  tip  may  be  seen  in  Figures  19  and  20. 

Time  histories  of  the  material  displacement  at  the  locations  indicated 
in  Figure  26  were  recorded  during  the  code  solution.  The  slippage  of  material 
at  three  points  along  the  crack,  as  given  by  the  distance  between  points  on 
opposite  sides  of  the  crack,  designated  as{W,X),  (R,S),  and  (M,  N)  in 
Figure  26,  is  shown  in  Figure  27.  The  material  on  the  right  side  of  the  crack 
is  moving  downward  and  to  the  left,  along  the  crack,  relative  to  the  material 
on  the  left  side.  The  extents  of  the  downward  (y-direction)  displacements  of 
these  points  and  the  crack  tip  (point  I)  as  a  function  of  time  are  shown  in 
Figure  28.  Tne  spatial  trajectory  of  the  point  pair  (W,X)  on  the  crack  surface 
during  the  time  span  of  the  solution  is  plotted  in  Figure  29.  Time-histories 
of  the  vertical  displacements  of  points  above  and  below  the  crack  tip  along 
four  vertical  cuts  (x  **  constant)  through  the  target  are  shown  in  Figures  30 
to  33.  The  downward  shift  of  material  persists  at  all  these  stations,  but  in 
smaller  amounts  as  the  distance  from  the  crack  increases. 

Stress  (o^)  -  time  profiles  at  points  H,  I,  J,  and  K,  along  the  hori¬ 
zontal  plane  through  the  crack  tip  (y  «  0),  are  shown  in  Figure  34.  Note  the 
increase  in  stress  over  that  of  the  loading  level  near  the  crack  tip.  Stress 
profiles  along  the  vertical  plane  through  the  crack  tip  (x  «  200  cm)  are  shown 
in  Figure  35.  The  points  above  the  crack  tip  show  a  peak  stress  of  5  kb  - 
corresponding  to  the  loading  or  incident  shock  ievel,  and  those  below  the 
tip  to  about  4  kb  -  the  transmitted  stress  level.  Stress-time  profiles  along 
a  vertical  plane  to  the  right  of  the  crack  tip,  at  x  ^  300  cm,  are  shown  in 
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Figure  $3.  Principal  Stress  Field,  Case 
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Figure  25.  Principal  Stress  Field,  Case 
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Figure  27.  Material  Slippage  at  Three  Points  Along  the  Crack,  Case  2. 
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Figure  29.  Spatial  Trajectory  cf  Two  Initially  Opposite  Points  on  the  Crack 
Surface ,  Case  2  . 
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Figure  31.  Vertical  Displacement  at  Points  Along  the  Vertical  Plane 
at  x  p~j200  cm.  Case  2. 
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Figure  32. 


Vertical  Displacement  at  Points  Along  the  Vertical  Plane 
at  x  %  250  cm ,  Case  2  . 
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Figure  33.  Vertical  Displacement  at  Points  Along  the  Vertical  Flane 
at  x  300  cm ,  Case  2 . 
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Figurt  3'’.  Stress  (ox)  Profiles  ;»t  Points  Along  the  Horizontal  Plane  at 
y  «  0,  Case  2. 
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Figure  36.  Here  the  stress  level  is  reduced  at  points  above  as  well  as  below 

the  crack  tip.  Time  histories  of  the  shear  stress  (o  )  at  points  along  the 

xy 

vertical  plane  through  the  crack  tip  are  shown  in  Figure  37. 

3.5  CASE  3  -  INTERACTION  OF  STRESS  WAVE  WITH  SINGLE, 

FINITE-LENGTH  CRACK,  WITH  CRACK  GROWTH 

For  this  case,  the  same  problem  as  in  Case  2  was  solved,  but  with 
the  provision  in  the  code  for  permitting  growth  of  the  crack  activated,  as 
described  previously. 

Representative  results  of  this  code  solution,  as  shown  by  the  particle 
velocity  fields  for  times  of  .5,  .6,  .7,  and  .9  msec,  are  given  in  Figures  5 
to  7  in  the  Summary,  Section  2.3.1,  and  in  Figure  30.  Associated  principal 
stress  fields  for  times  of  .5  and  .6  msec  are  shown  in  Figures  39  and  40. 
Propagation  of  the  crack  does  occur  for  this  loading;  the  extent  of  the  crack 
at  any  time  is  indicated  by  the  circled  lattice  points.  For  this  loading  function, 
the  fracture  criterion  used  was  met  at  relatively  late  times  after  the  shock  front 
has  passed,  so  that  the  effect  of  crack  growth  on  the  transmitted  wave  is 
probably  not  large.  These  and  other  effects,  such  as  material  slippage,  will 
be  examined  through  comparisons  of  recorded  time  histories  for  Case  2  and 
Case  3  in  forthcoming  analyses. 

4.  DYNAMIC  GRIFFITH  CRITERION  AND  CRACK  PROPAGATION 

To  improve  the  physical  significance  of  the  numerical  solutions,  the 
conclusion  was  reached  that  a  more  realistic,  dynamic  criterion  for  predicting 
and  following  the  course  of  crack  propagation  in  flawed,  jointed,  brittle  media 
should  be  incorporated  in  the  numerical  method.  A  necessary  step  in  making 
such  a  formulation  change  is  to  check  the  new  code  through  comparisons  with 
analytical  results  from  some  model  problems,  as  will  be  discussed  in  Section  5. 

The  modifications  needed  to  follow  crack  propagation  are  threefold. 
These  are  the  logic  associated  with  the  crack  geometry  and  crack  propagation, 
the  velocity  threshold  criterion,  chosen  here  to  be  the  Griffith  criterion,  and 
the  velocity  orientation  criterion,  chosen  to  be  the  direction  transverse  to  the 
maximum  principal  stress  vector.  Programming  of  tills  formulation  is 
currently  underway. 
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Figure  36.  Here  the  stress  level  is  reduced  at  points  above  as  well  as  below 

the  crack  tip.  Time  histories  of  the  shear  stress  (o  )  at  points  along  the 

xy 

vertical  plane  through  the  crack  tip  are  shown  in  Figure  37. 

3.5  CASE  3  -  INTERACTION  OF  STRESS  WAVE  WITH  SINGLE, 

FINITE-LENGTH  CRACK,  WITH  CRACK  GROWTH 

For  this  case,  the  same  problem  as  in  Case  2  was  solved,  but  with 
the  provision  in  the  code  for  permitting  growth  of  the  crack  activated,  as 
described  previously. 

Representative  results  of  this  code  solution,  as  shown  by  the  particle 
velocity  fields  for  times  of  .5,  .6,  . 7,  and  .9  msec,  are  given  in  Figures  5 
to  7  in  the  Summary,  Section  2.3.1,  and  in  Figure  38.  Associated  principal 
stress  fields  for  times  of  .5  and  .6  msec  are  shown  in  Figures  39  and  40. 
Propagation  of  the  crack  does  occur  for  this  loading;  the  extent  of  the  crack 
at  any  time  is  indicated  by  the  circled  lattice  points.  For  this  loading  function, 
the  fracture  criterion  used  was  met  at  relatively  late  times  after  the  shock  front 
has  passed,  so  that  the  effect  of  crack  growth  on  the  transmitted  wave  is 
probably  not  large.  These  and  other  effects,  such  as  material  slippage,  will 
be  examined  through  comparisons  of  recorded  time  histories  for  Case  2  and 
Case  3  in  forthcoming  analyses. 

4.  DYNAMIC  GRIFFITH  CRITERION  AND  CRACK  PROPAGATION 

To  improve  the  physical  significance  of  the  numerical  solutions,  the 
conclusion  was  reached  that  a  more  realistic,  dynamic  criterion  for  predicting 
and  following  the  course  of  crack  propagation  in  flawed,  jointed,  brittle  media 
should  be  incorporated  in  the  numerical  method.  A  necessary  step  in  making 
such  a  formulation  change  is  to  check  the  new  code  through  comparisons  with 
analytical  results  from  some  model  problems,  as  will  be  discussed  in  Section  5. 

The  modifications  needed  to  follow'  crack  propagation  are  threefold. 
These  are  the  logic  associated  with  the  crack  geometry  and  crack  propagation, 
the  velocity  threshold  criterion,  chosen  here  to  be  the  Griffith  criterion,  and 
the  velocity  orientation  criterion,  chosen  to  be  the  direction  transverse  to  the 
maximum  principal  stress  vector,  programming  of  this  formulation  is 
currently  underway. 
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Figure  40.  Principal  Stress  Field,  Ca.i 


It  was  suggested  by  Griffith  that  a  useful  way  of  accounting  for 
controlled  crack  growth  in  an  elastic  body  is  to  balance  the  rate  of  work 
(power  input)  to  the  body  by  the  rate  of  uptake  of  strain  energy,  kinetic 


energy,  and  surface  energy,  i.e., 

P  =  Uj  +  V2  +  U3  (4) 

where 

P  =  S  lij  vj  d  ai  (5) 

Ui=/WdV0  (6) 

U2=%Jvk''kdVo  <7> 

U3  M  r  $  nk  dak  (8) 


W  is  the  strain  energy  density,  T  is  the  specific  surface  energy,  a  material 
constant,  and  the  dots  denote  material  differentiation.  The  material  deriva¬ 
tives  consist  of  two  components,  i.e., 


f  =  * 

1  at 


X,  6 


af_ 

as 


L  « 

x,  t 


(9) 


where  6  is  a  measure  of  the  crack  length.  When  the  crack  is  not  moving  f 
is  the  ordinary  material  derivative. 

In  order  to  accommodate  the  crack  motion  it  is  necessary  to  follow 
this  procedure.  At  each  point  in  time,  a  test  iteration  is  carried  out  in 
which  the  crack  is  assumed  to  run  a  small  distance.  New  (Uj,  U3 ,  U^Jare 
calculated  for  the  new  6,  and  Eqn.  (4)  is  entered  and  checked.  If  the  right 
hand  side  is  less  than  the  left  hand  side,  the  crack  either  has  not  yet  started 
to  run  or  the  assumed  increase  in  6  was  too  large.  We  propose  to  choose  a 
small  enough  value  of  6,  such  that  the  error  incurred  by  always  assuming  the 
crack  has  not  run  when  the  inequality  occurs,  is  a  tolerable  error.  When 
equality  occurs,  the  crack  is  advanced  another  small  increment,  and  the 
program  continues. 


56 


I 
i : 
€ 


s 

i 


“■  *y-  > tv ‘/'Mm  •-  -•  ^ *» '< £ynzp,$j~ btti&i**'! *' ? 1  ,>y>'i.«4v- 


At  each  stage  the  program  will  calculate  the  direction  of  the  principal 
maximum  stress  vector  and  the  crack  will  move  in  a  direction  transverse  to 
that.  The  crack  geometry  will  be  coded  entirely  like  a  free  surface  with  the 
same  boundary  type  characterization  that  is  currently  being  used. 

5.  ANALYTICAL  COMPARISON  PROBLEMS 


The  only  elasto-dynamic  solutions  that  are  currently  available  for 

7  8  9 

an  accelerating  crack  are  those  of  Kostrov  ,  Eshelby  ,  and  Achenbach  .  The 
crack  growth  is  produced  by  anti-plane  shear.  It  was  thus  decided  to  modify 
the  SHEP  code  to  accommodate  anti-plane  shear  or  out-of-plane  displacement 
in  a  manner  that  would  retain  the  feature  that  the  out-of-plane  displacement 
remains  independent  of  the  z-coordinate  (or  S-coordinate  for  the  axisymmetric 
case).  Thus,  the  additions  maintain  the  two-dimensional  character  of  SHEP 
because  the  out-of-plane  motion  merely  superimposes  on  and  does  not  couple 
with  the  in-plane  motion,  assuming  that  we  are  confined  to  the  linear  elastic 
region. 

The  additional  equations  that  have  been  incorporated  into  the  SHEP 
code  to  account  for  out-of-plane  motion  are  shown  below  in  differential  form. 
The  difference  forms  follow  the  same  format  as  used  by  the  rest  of  SHEP. 

We  have 


dr 

__£z 

dy 


r 

yz 


=  p  W 


(10) 


E  -  (P+  q)  V  +  V  (s^  ex+  syey  +  s0e0  4  Txy  exy  +  Txz  Scz  +  ryz  ®yz^ 

(ID 


xz 

=  G  e  xz 

(12) 

yz 

=  G  e 

yz 

(13) 

xz 

n  dw 
=  G  dT 

(14) 

yz 

r  dw 

=  G  sF 

(15) 
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All  the  rest  of  the  equations  retain  the  same  form  (the  form  of  the  Jauman-Oldroyd 
derivatives  is  given  in  the  Appendix) . 

A  test  problem  consisting  of  simple  shear  motion  of  a  slab  has  already 
been  successfully  run  with  the  modified  code,  as  discussed  in  the  Appendix. 

The  next  car,e  that  will  be  run  involves  the  loading  of  a  stationary  crack  in 
anti-plane  shear. 
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APPENDIX 

CORRECTION  TO  THE  ROTATIONAL  TERM  IN  THE  STRESS  CALCULATIONS 

During  the  examination  of  the  formulation  of  the  code  preparatory  to 

making  the  modifications  discussed  in  Sections  4  and  5  of  the  report ,  an  error 

in  the  sign  of  the  rotational  correction  term  in  the  stress  calculations,  as 

A1 

originally  proposed  by  Wilkins  ,  was  discovered.  In  addition,  from  the 

A2 

standpoint  of  material  objectivity  ,  the  form  of  the  rotational  term  is  si  on 
to  be  written  in  only  an  approximate  form.  These  statements  are  documented 
below.  The  sign  was  corrected  in  our  code  and  a  check  of  the  correction 
was  made  by  running  a  model  problem.  The  results  confirmed  the  sign  correc¬ 
tion  and,  in  addition,  verified  the  accuracy  of  the  SHEP  code. 


In  SHEP,  the  constitutive  equations  are  correctly*  subsumed  in  the 


formA3: 

*  su  -  6ij  p  ® 

(Al) 

% 

=  2G  (By  -  j  ®  6ij  ) 

(A2) 

0 

=  In] 

(A3) 

J 

dV 
=  dVo 

(A4) 

where 

liJ 

is  the  cartesian  Couchy  stress  tensor 

sij 

is  the  stress  deviator  tensor 

P 

is  the  pressure 

J  is  the  local  volume  ratio,  equal  to  the  square 
root  of  the  third  stretch  invariant. 


♦The  term  "correctly"  implies  that  Eqn.  (A2)  satisfies  the  principle  of 
material  objectivity. 
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Preceding  page  blank 


G  is  the  shear  modulus 

0  is  the  dilatation 

• 

is  the  Cartesian  strain  rate  tensor 

V  is  the  local  instantaneous  specific  volume 

and  in  (A4) , 

is  the  particle-differentiator. 

From  Eqn.  (A3)  and  (A4),  it  follows  that 

'-H 

II 

<|<* 

(AS) 

• 

Furthermore  .  n  Eqn.  (A2),  the  dot  over  e^  denotes  the  material  derivative 
while  the  hat  over  s^  denotes  the  Jaumann-Oldroyd  derivative.  The  form 
of  the  latter  is  given  by: 

=ik  -  Jik  -  vi,j ‘jk +  Vj,k 

• 

(A6) 

where 

• 

vi.i  *  e«  +  “ij 

(A7) 

or 

e 

2  -  ij 

(A8) 

where 

OJjj  is  the  spin  tensor. 

The  last  two  terms  on  the  right  hand  side  of  Eqn.  (A6)  correspond  to 
Wilkins'  6jj;  a  simple  calculation  reveals  that  the  entire  form  of  Eqn.  (B— 14) 
on  page  78  of  Reference  A1  is  formally  incorrect  but  approximately  appropriate 
(apart  from  a  sign  error). 

The  sign  error  apparently  arose  in  the  following  way:  Wilkins  intro¬ 
duces  an  angular  velocity  vector,  which  he  denotes  by  [sin  w  ].  The  angular 
velocity  vector  is  correctly  taken  as  the  curb  of  the  velocity  vector  which, 
according  to  the  usual  right-handed  correction,  corresponds  to  a  counter- 
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clockwise  rotation.  Wilkins  then  uses  the  rotation  angle  to  to  transform  the 

stress  components  to  new  rotated  components.  The  transformation  equations 

A3 

he  adopts  are  taken  from  Timoshenko  .  Presumably  Wilkins  changed  the 
signs  of  the  terms  that  are  odd  in  to,  because  Timoshenko  shows  a  diagram 
in  which  the  rotation  is  clockwise.  What  one  must  notice ,  however,  is  that 
Timoshenko's  picture  is  based  on  a  left-handed  system  of  axes.  Thus  Timo¬ 
shenko's  equation  should  have  been  subsumed  without  the  sign  change. 

To  further  check  our  logic,  we  ran  the  following  simple  problem.  A 
slab  is  set  into  simple  shear  motion,  of: 


Nonlinear  (Poynting)  effects  are  associated  with  the  lengthening  of  the 
originally  upright  fibers.  The  correct  rotation  terms  (with  the  right  sign)  pre¬ 
dict  a  tensile  stress  in  the  lengthened  fibers,  in  agreement  with  the  analytical 
solution.  The  original  program,  with  the  incorrect  sign,  predicts  compression. 
Keep  in  mind  that  these  are  second  order  stresses  and  thus  the  effect  of  this 
error  on  the  overall  computation  becomes  important  only  when  the  first  order 
stresses  are  comparable  to  the  shear  modulus  of  the  material.  The  reason  is 
that  the  second-order  stresses  go  as  the  square  of  the  first-order  stresses, 
e.g. ,  in  the  shear  problem: 


(A9) 


In  the  process  of  making  these  changes  in  SHEP,  we  checked  out  a 
model  problem,  namely  the  simple  shear  at  constant  velocities  of  an  infinitely 
long  slab  (a  1-D  problem).  The  results,  in  terms  of  the  stress  at  the  moving 
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surface  and  the  stress  at  the  fixed  surface  are  plotted  in  Figures  Al  and  A2, 
respectively.  The  analytical  solution  is  just  a  sum  of  Heaviside  step  functions. 
The  agreement  witr.  the  analytical  solution  is  attested  to  by  the  fact  that  the 
computed  jumps  fall  directly  (within  1%)  on  the  grid  line  (both  coordinate  axes 
are  non-dimensionalized) .  The  general  solution  looks  like 


where 


T  ® 

-  —3—  =  L  L  H(s  -  y  -  2nh]  +  H  (s  +  y  -  2h  -  2nh)  ]  (A10) 

Mo  n  =  0 

v 

MQ  (the  Mach  number)  =  — - 
s  =  ct 


and  h  is  the  height  of  the  slab. 
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Figure  Al.  Shear  Stress  at  the  Moving  Surface 


rw*  ■  i 
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